function [F,p] = ftest(r,fg,npoints,varargin)
% Input:  r  = list of squared sums of residuals*
%         fg = list of degrees of freedom (fit parameter)*
%         npoints = number of data points
%         alpha (optional) = value of significance (default = 0.001)
% Output: p = list of p-Values
% * Start with simpliest model
if nargin == 3
    alpha = 0.001;
elseif nargin == 4
    alpha = varargin{1};
end
d(1:length(r)) = npoints-fg(1:length(r));
for i=2:length(r)
if(fg(i-1)>fg(i)) disp('Start with simpliest model');break; end
F(i-1) = ( (r(i-1)-r(i) )/ (d(i-1)-d(i) ) ) / ( r(i)/d(i) );
p(i-1) = 1-fcdf(F(i-1),d(i-1)-d(i),d(i));
if p(i-1)<=alpha
    fprintf('Fit %i is better than Fit %i (p = %d, F = %d)\n',i,i-1,p(i-1), F(i-1))
elseif p(i-1)>=alpha 
     fprintf('Fit %i NOT better than Fit %i (p = %d, F = %d)\n',i,i-1,p(i-1), F(i-1))
end
end